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ABSTRACT 



Methods used in the radial-velocity program of short-period binary systems at the 
David Dunlap Observatory are described with particular stress on the Broadening Func- 
tion (BF) formalism. This formalism has permitted determination of radial velocities 
from complex spectra of multiple-component systems with component stars showing 
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' very different degree of rotational line broadening. The statistics of random errors of 

orbital parameters is discussed on the basis of the available orbital solutions presented 
in the six previous papers of the series, each with ten orbits. The difficult matter of sys- 
tematic uncertainties in orbital parameters is illustrated for one typical case of GM Dra 
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^) . from the most recent Paper VI. 

Qh' Subject headings: stars: close binaries - stars: eclipsing binaries - stars: variable stars 

o 

1. INTRODUCTION 

^ ' This paper should be considered as a companion and supplement to the previous papers of our 

' series of radial velocity studies of close binary stars: Lu &: Rucinski (1999, Paper I), Rucinski k, 

Lu (1999, Paper II), Rucinski et al. (2000, Paper III), Lu et al. (2001, Paper IV), Rucinski et al. 
(2001, Paper V), Rucinski et al. (2002, Paper VI). 

The current program of radial velocity observations of close binary systems with periods shorter 
than one day is approximately at its half-way point. Our methods have been evolving slightly during 
the execution of the 60 radial velocity orbits presented in the six papers of the series, but appear to 
have stabilized now, warranting a more detailed documentation of the essential steps in our analysis 
and data reductions. We summarize these methods and give an overview of the uncertainties so 
that the results described in the previous and the planned future papers of the series could be 
better evaluated by readers. The discussion is limited strictly to methodological aspects and does 
not include any astrophysical results which will be discussed after the program is concluded. 



^Based on the data obtained at the David Dunlap Observatory, University of Toronto. 
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2. INSTRUMENTATION AND OBSERVATIONS 

We observe radial velocities of close binary stars with the 1.88m telescope of the David Dunlap 
Observatory using its medium-resolution spectrograph in the Cassegrain focus. The angular scale in 
the telescope focus is 6 arcsec per mm. We use one of the two spectrograph slits, 300 /xm or 250 /xm 
in width, both fixed in the E-W orientation and both 10 mm long. The angular widths of 1.8 and 
1.5 arcsec approximately match the median seeing at the DDO of 1.7 arcsec. Since we started with 
the shortest-period binaries showing the strongest rotational line broadening, most observations 
have been made with the 300 /xm slit. The scale reduction of the collimator-camera combination 
is 4 times resulting in the slit image of 75 /xm or 62 for either of the slits. Our light detector 
is currently a thick, front-illuminated CCD chip of 1024 by 1024 pixels, 19 /xm square. Thus, the 
slit images have the total widths of 3.9 or 3.3 pixels while the FWHM widths arc about 2.6 and 
2.2 pixels for the respective slits. To lower the influence of the readout noise, the two-dimensional 
CCD images are on-chip binned 4-times in the direction perpendicular to the dispersion direction. 

Most of the spectral data have been obtained using the 1800 lines/mm diffraction grating, with 
the spectral window centered on the magnesium triplet Mg I "b" at 5184 A. For solar-type stars, 
this region is very rich in spectral lines, which is an essential consideration for our method of radial- 
velocity measurements - through broadening functions - to succeed. The main-sequence stars of 
spectral types of middle-A to middle-K are practically the only stars found in close binaries with 
orbital periods shorter than one day. At 5184 A, the spectrograph delivers 0.204 A/(19 fim pixel) 
or about 11.8 km s^^ per pixel. As it is well known, when cross-correlation or similar techniques 
are used, narrow, properly-sampled, symmetric spectral features can be usually measured to better 
than about 1/10 part of the pixel size, with the accuracy growing in relation to the total length 
of the spectrum. In our case, the spectrum has the length of 208 A so that we can rather easily 
determine velocities of sharp-line stars with accuracy of about 1 km s~^, as has been verified by 
many observational programs at the DDO (see the end of this section). The accuracy for broad-lined 
spectra of binary components is obviously lower and depends on a combination of many factors. 
We discuss the random errors in Section 8, while systematic uncertainties specific to close binary 
stars are discussed in Section 9. 

The spectrograph is known to show some flexure limiting exposures to about 30 minutes. This 
has not been a real limitation in the observations because our program stars have short periods 
and in fact require exposures to be no longer than 15 - 20 minutes to prevent the radial-velocity 
smearing. We normally take comparison-lamp (FeAr) spectra before and after each exposure, but 
sometimes, for the shortest-period binaries, we take two or three stellar exposures for a set of 
bracketing comparison spectra spaced by no more than half an hour. The flat field lamp is an 
internal one in the spectrograph, but we occasionally take also sky flat-field spectra. 

We observe typically 3 to 5 radial- velocity standard stars per night. These stars are selected 
to have spectra of similar spectral types to our program stars, to serve later as templates in our 
technique of radial- velocity measurements through broadening functions (see Sections 5-6). An 
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inter-comparison of radial velocities of standard stars gives an estimate of random errors at the 
level of about 1.0 - 1.2 km s'^. This agrees with the results for Cepheids observed at the DDO by 
Sugars &: Evans (1994) and Evans (2000) where the errors were estimated at 1.0 — 1.3 km s~^. A 
fraction in these errors may come from our continuing use of the lAU standard-star list as published 
in the Nautical Almanac 1995. As explained in Stefanik et al. (1999), the lAU list contained at 
that time a few stars which are unsuitable as radial velocity standards. Since we used many 
different standard stars from the lAU list, these systematic errors averaged out to some degree and 
manifested themselves mostly as random errors. We now use exclusively the list of Stefanik et al. 
(1999), but the results of our program may be affected by the uncertainties in the old standard 
velocity data at a level of 0.2 - 0.5 km s~^. 



3. INITIAL ANALYSIS OF SPECTRA 

The reductions consist of several stages. Stage #1 consists of a transformation from two- 
dimensional images to one-dimensional, wavelength-calibrated, rectified spectra. All the steps, 
starting with de-biasing and flat-fielding, utilize the standard techniques within IRAF^. We make 

sure to use a consistent set of low-order polynomials for the dispersion relation and use the IRAF 
rejection algorithm for rectification of the spectra; both steps are facilitated by the short length of 
our spectra within which the dispersion and the CCD sensitivity vary only slightly and in a smooth 
way. 

Figure 1 shows a typical spectrum for our program of the binary KR Com, one of the stars 
presented in the most recent paper of the series (Rucinski et al. 2002, Paper VI). It is a typical, 
yet difficult case, in the sense that we have been frequently dealing with rather complex, multi- 
component spectra even among relatively bright stars (7th magnitude in this case). Some systems of 
our program have been previously known binaries, too difficult to handle using traditional (including 
cross-correlation) methods, but many were recently discovered as photometric variables. 

The spectra of KR Com are dominated by the third, slowly rotating component which - 
although the fainter one in the visual system - dominates the spectral appearance and produces 
sharp, easily identifiable spectral lines. The triplicity of KR Com went apparently unnoticed for 
so long mostly because the spectra, superficially, look like those of a single, slowly-rotating star 
and - paradoxically - the spectrum of the brighter binary component is not normally visible. The 
presence of the binary, which produces the broad spectral signature, manifests itself spectrally only 
through merging of more common, weaker lines and the general lowering of the continuum and is 
difficult to notice in low S/N spectra. The low-level photometric variability of the whole system is 
due to the contact binary which is the brighter component in the system, but the variability signal 



^IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of 
Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. 
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is sufficiently "diluted" in the combined light of the system that it took the high quality of the 
Hipparcos photometry to discover it. 

Spectra like that shown in the left panel of Figure 1 are not analyzed directly, but are subject to 
the broadening-function extraction process, which is followed by measurements of radial velocities. 

4. WHY WE DO NOT USE THE CROSS CORRELATION FUNCTION 

Step 7^2 of the analysis is determination of the broadening function (BF). The BF approach 
was described before in Rucinski (1992, 1999) and is discussed more extensively in Section 5. In 
essence, it consists of a linear, least-squares determination of the broadening convolution kernel 
from rotationally- and orbitally-broadened spectra, utiUzing spectra of sharp-line, slowly rotating 
radial-velocity standards. We do not use the popular cross-correlation function (CCF) technique 
because it appears to give inferior and biased results for close binary systems. We now try to 
explain this rather strong statement. 

1. The CCF combines the broadening of the program spectrum with that of the template, with 
the resulting loss of resolution, while the BF approach attempts to remove the common 
broadening contributions. Only if the template spectrum were a series of delta-functions, 
would the results be the same. 

2. The definition of the baseline in the CCF is usually difficult and may lead to problems when 
relative luminosities of components are determined. 

3. Outside of the main peak which is used for radial velocity determination, the CCF always 
shows a fringing pattern which may affect the strength and intensity of secondary correlation 
peaks for multiple systems. For very close binaries, the secondary fringes frequently produce 
the "peak-pulling" effect of the systematically smaller radial- velocity amplitudes, but a more 
complex interaction is entirely possible. 

4. The shape of the CCF beyond the main correlation peak depends on the shape of the stellar 
spectrum. For the same star, observations in different parts of the spectrum define different 
CCF's. This problem is rarely recognized and is particularly severe for sparse spectra, when 
the CCF is analyzed over a wide range of correlation lags. 

The problems listed above are illustrated in Figure 1 for the case of the triple system KR Com. 
In the right panel, we show a comparison of the BF with the CCF for the same spectra. While the 
BF very clearly shows all three components in the system, it would be very difficult separate the 
three signatures using the CCF. The superior resolution offered by the BF approach has permitted 
us to analyze spectra with very strong rotational broadening, combined with situations of three or 
more sets of blended lines in triple and quadruple systems, with component stars showing different 
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amounts of rotational broadening. Such systems have been frequently abandoned in the past 
because of insurmountable difficulties with separating and measuring radial velocities of individual 
components. 

The problems of the baseline location and fringing in the CCF, as well as of the dependence 
on the shape of the stellar spectrum are illustrated in Figure 2. This figure contains a result of 
the following data-processing experiment. A high quality, but somewhat sparse stellar spectrum 
(left panel, sampled at equal velocity steps of 0.88 km s"-*^) was convolved with the single-star 
rotational-broadening pattern with Vsini = 88 km and then subject to the CCF determination. 
No noise was added, so the CCF is basically perfectly determined and can be used for a direct 
comparison with the assumed broadening function. The right panel of the figure compares the 
assumed broadening profile (dotted line), which is the same as the BF, with the CCF (full line). 
The strong fringes in the CCF are very well visible. In this particular case, the unusual strength 
the positive fringes results from the low density of spectral features in the original spectrum and 
illustrates the dependence of the CCF on the stellar spectrum. While the BF formalism is insensitive 
to the density and distribution of the spectral lines, the CCF - beyond the main peak - does depend 
on the spectral region. Thus, sparse spectra will lead to less-well defined broadening functions 
with larger random errors (simply because of the lower information content), while the CCF will 
additionally show systematic differences in the fringing pattern outside the main peak. 

The negative fringes which are always present in the CCF can produce a quasi-bascline around 
the main correlation peak at a very different level than expected. In the case shown in Figure 2, 
the local baseline is located about —0.1 below the originally assumed broadening profile. If only 
a small part of the CCF were used, this is where the baseline would normally be located. Since a 
CCF would rarely be used for anything else but a radial-velocity determination from the correlation 
peak, the exact location of the baseline may seem immaterial. However, when a secondary star is 
added to the picture, with the similar rotational broadening and a velocity separation comparable 
to the rotational broadening, as is the case for short-period binary stars, then the secondary peak 
in the CCF will definitely interact with the primary-star fringing pattern; there will be also the 
reverse interaction of the secondary pattern with the primary peak. For single objects, the fringes 
are basically irrelevant, so that very close or identical radial velocities are determined using both 
techniques. Problems occur when multiple components are present in the spectra and are partic- 
ularly severe in the CCF when broadening of the lines is of the same scale as the line splitting, 
exactly the situation we face in our program. 

We recognize that a method based on the two-dimensional cross-correlation function called 
TODCOR (Zucker & Mazeh 1994; Zucker et al. 1995) has been developed and successfully applied 
to several multi-lined stellar systems showing sharp spectral lines. We did not attempt to use this 
technique mostly because we feel more comfortable with a tool developed by ourselves, but also 
because (1) TODCOR is designed for sharp-lined spectra and has not been demonstrated to work 
for the very broad lines of contact binaries, and (2) we frequently deal with mixed very broad 
and narrow spectral signatures which would require extension of the TODCOR capabilities even 



-6- 



further. We note that the non-lineax nature of the cross-correlation compUcates derivation of relative 

luminosities of components and requires a complex calibration while our linear approach gives 
directly the relative luminosities through integration of the individual features in the broadening 
functions. This is particularly convenient for systems with components showing very different 
degrees of rotational broadening. 

5. BROADENING FUNCTIONS 

We define the broadening function^ as a function that transforms a sharp-line spectrum of a 
standard star into a broadened spectrum of a binary or, for that matter, of any other star showing 
geometrical, Doppler-effect line broadening. This way we not only determine the broadening- 
function shape, but also automatically relate absolute velocities of program stars to radial- velocity 
standards used as templates, a common advantage with the CCF approach. We do not use model 
spectra, e.g. through representation of spectral lines by delta functions. While the broadening 
functions determined that way would be cleaner and much better defined than those utilizing 
standard-star templates, the advantage of the automatic relative radial-velocity calibration would 
be lost. 

We perform all radial velocity determinations in the geocentric system and only later transform 
the results to the barycentric (heliocentric with planetary corrections) system. Thus, we start with 

a raw template spectrum St, with its wavelength scale in Wt, and a raw program spectrum Sp, 
with its wavelength scale in Wp. Both St and Sp arc rectified and normalized to unity. To diminish 
importance of the one-to-zero discontinuities at the ends of the spectra, we invert them so that the 
absorption lines are represented by positive spikes: S^ = 1 — St and Sp = 1 — Sp. 

The spectra must be of similar spectral type. We normally use the templates with spectra 
within one spectral type. However, we have found that F-type templates will work reasonably 
well for radial velocity determinations between middle A-types to early K-type stars; however, the 
relative luminosity estimates from the individual peaks will then be wrong. 

The spectra must be first re-sampled into equal steps in velocity. In our case the velocity step 
is typically Av = 11.8 km s"^. An auxiliary vector of wavelengths is now created with elements: 
W = Wq * (1 + r)*, where i = 0, . . . ,n — 1 is the index in the new vector and r = Av/c, where c is 
the velocity of light. The origin of this vector, Wq, is selected to fall just above both origins of Wt 
and Wp for a meaningful interpolation of both spectra into the new wavelength scale. The length 
of W is in our case usually selected to be n = 1000 — 1020 spectral elements. The spectra S'^ and 
Sp are linearly interpolated using W, by treating Wt and Wp as the respective abscissas, to create 
the spectra used in the BF derivation: for the template, T, and for the program star, P. After this 



full description of the concept of the broadening functions with examples and detailed programming suggestions 
is available at the WWW site http://www.astro.utoronto.ca/~rucinski. 
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is accomplished, the three wavelength vectors W, Wt and Wp are no longer needed because the 

program and the template spectra are now in the same (geocentric) velocity system. We can think 
about them as functions P(n) and T(n) with the same velocity axis or vectors P and T over the 
same range of indices. 

The convolution operation which maps a sharp- line spectrum T into a broad and/or binary-star 
spectrum P, 

P(A') = J B{\' - A) r(A) d\ (1) 
can be written as an array operation, 

P = DB (2) 

in which the rectangular array D is created from the vector T by placing it as columns of D 
after shifting it downward by one index for each successive column (sec below or, for further 
details, consult Rucinski (1992, 1999)). The broadening function is represented by a vector of the 
unknowns in the solution, B. The array D has the short dimension m and the long dimension 
n — m + l; it accomplishes the mapping of T ^ P. We normally use the odd number for the size of 
the broadening function, m, to have it centered at the pixel symmetrically distant from both ends. 
Also, for proper handling of the ends, m' = integer (m/2) points are removed from both ends of P. 

The convolution operation equivalent to Eq. (2), which is used in the least-squares determina- 
tion of B, can be written as a system of over-determined linear equations: 

m~l 

Pi = Ti^jn-j Bj with i = m', . . . ,n — m' — 1 (3) 

i=o 

The number of equations should be several times larger than the number of unknowns, n — m + l > 
rn. In our program, wc normally use n = 1000 — 1020 and m = 121. The size of the broadening 
function, m, translates into the relative velocity range (program minus template) of ±708 km s^^ 
insuring a good definition of the BF itself and of the fiat baseline around it. The actual size of 
the broadening function is a matter of choice; sometimes we repeat the BF determination with a 
smaller m for binary systems with moderate line splitting when a wide window of over 1400 km s^^ 
is not needed. The point is to use as short a BF as possible because the quality of the determination 
(over-determinacy) increases in relation to how many times the spectra are longer than the BF. 

Solving of the broadening function Bj is accomplished by least-squares. We are strong advo- 
cates of the Singular Value Decomposition (SVD) technique which is particularly useful in eliminat- 
ing those parts of the spectra that carry no information (the inter- line continuum) , but create linear 
dependencies. The approach involving rejection of small singular values is the best for restoration 
of the shape of the BF for its subsequent modeling. However, with the radial velocities in mind, we 
do not in fact eliminate any singular values. In this respect, we have made some departure from 
the original philosophy, but this departure has a reason; if some basis functions are eliminated, 
there exists a possibility that the spectral features may acquire asymmetries through an unwanted 
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conspiracy of the basis functions which remain in the definition of the BF. By retaining ah singular 

vahies, we treat each clement of Bj (Eqs. (2) - (3)) as a totally independent variable, not related 
in any way to its neighbors. Thus, any least-squares technique can be used at this stage, although 
we continue to use the SVD as easier to use and more transparent for the matrix inversion. 

The details of the SVD approach to solve the array equation Eq. (2) or its equivalent Eq. (3) 
for the BF vector B are described in Rucinski (1992) and the programming examples are given in 
Rucinski (1999). Even without elimination of any singular values in the SVD solution, this approach 
has an advantage that one decomposition of the template-spectrum array, D = U W V""", can serve 
to determine B from several program spectra through the inverted relation, B = VW~^ (U'^P). 
For an excellent exposition of the SVD technique stressing its beneficial properties, see Press et al. 
(1992). 

Irrespective of which method of the least-square solution is used, the resulting broadening 
functions are always very noisy and cannot be used for radial-velocity measurements. The reason 
for the excessive noise is that each element of the solution Bj is unrelated to its neighbors and 
is treated as a separate unknown. We know, however, that our spectral resolution is controlled 
by the spectrograph slit which introduces coupling between neighboring points of the BF. In our 
case, the intrinsic smoothing introduced by one of the entrance slits is characterized by the FWHM 
of about 2.6 or 2.2 pixels. It is therefore reasonable to apply some smoothing to the noisy BF's. 
Superficially, this step does the same to the final shape of the BF as smoothing through rejection 
of noise absorbed by high-order singular values in the SVD technique; however, this operation is 
strictly local whereas removal of some singular values may introduce non-local effects. Usually, 
we smooth the broadening functions by convolving them with a Gaussian with a = 1.5 pixels 
(FWHM=3.53 pixels); for poor spectra of very faint stars we are sometimes forced to use a = 2.0 
(FWHM=4.71 pixels). Such smoothing is slightly stronger than its instrumental counterpart by 
the spectrograph slit, but is nevertheless very small, when compared with widths of lines in binary 
stars with periods shorter than one day. 



6. RADIAL VELOCITY MEASUREMENTS 

Step #3 is the radial velocity determination from the broadening functions. We determine 
the radial velocities of each binary component in the geocentric system, relative to the template 
star, Vi- Following that, the relative velocities are transformed to the solar system barycenter with 
Vi = Vi + {HCp — HCt) + Vt, where HC arc the barycentric (sometimes called "heliocentric") 
velocity corrections resulting from the orbital Earth motion for the program and template spectra 
and Vt is the barycentric velocity of the template star. 

The broadening function B{v) determined in the previous step is defined at points separated 
by equal steps in relative geocentric velocity, in our case normally At; = 11.8 km s^^, spanning 
the velocity range —708 <V< -1-708 km s~^. The velocities of stellar components are determined 
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by simultaneous fitting of several Gaussian curves to that many spectral features as are seen in 
the BF. Thus, it would be a 4-parameter Gaussian fit for a single star (baseline oq, strength ai, 
position 02, width 03), a 7-parameter fit for a binary or a 10-parameter fit for a triple system, etc., 
as in: 

n 

B{v) ~ 00 + ^ au exp {{v - a2i)/asif^ (4) 

i=l 

where n represents the number of stellar components in the system. We found that least-squares 
Gaussian fits for single stars are usually stable, while those involving more components (binary, 
triple and higher multiplicity systems) are numerically unstable forcing us to fix or manually adjust 
the width parameters, a^i. 

In triple systems which we encountered so far, the most typical combination has been a broad- 
lined close binary accompanied by a sharp-lined, slowly-rotating star. In such situations, we first 
leave the width and position of the third, sharp component floating in order to determine best 
possible parameters for its subsequent subtraction from the BF. For the example shown in Figure 1, 
the Gaussian widths for the binary components were assumed at 031 = 110 km s~^ and 032 = 70 
km while the width 033 was determined at 24.74 km s~^. We found that situations similar to that 
shown in the figure require a careful removal of the third-component signature. In order to define 
the BF for the close binary the best way possible, we cannot remove the averaged signature of the 
third star from many spectra and must subtract it as it is defined for the same observation. There 
may be many reasons why subtraction of the averaged third-star peak leaves too large residuals, 
including small changes in the effective resolution, imperfections in the geocentric to barycentric 
transformations or instabilities in the spectrograph. Obviously, this approach reduces the accuracy 
of the third-star velocities, but our goal has been to determine the best velocities for the close 
binary so that we accept this limitation. The BF for the binary is usually very well defined, see for 
example Figure 4 in Paper IV for HT Vir. 

The random radial- velocity errors for binaries in occurring in triple systems are only slightly 
larger than for the isolated binaries, typically by less than 1 km s~^ in the errors of Vq, Ki and K2, 
which simply reflects more degrees of freedom in the problem (see the discussion in Section 8 and in 
Figure 4). Much more difficult to characterize are systematic uncertainties. One manifestation of 
such uncertainties is the occasional presence of an undesirable "cross-talk" in the three-feature fits in 
the sense that the third-star velocities sometimes correlate with the binary phase. We always check 
the third-star velocities for the dependence on the binary phase and sometimes find a correlation. 
(For a rather extreme case, see the discussion of II UMa in Paper VI). It is usually quite difficult 
to find reasons for the cross-talk and each case seems to be unique. Faintness of the star and/or 
poor spectra certainly magnify the problem which depends on such factors as the location of the 
third peak in the BF relative to the peaks for the binary system stars (i.e. with which component 
the third peak mergers most of the time) or the overall degree of the line splitting for the binary 
system (which depends on the orbital inclination). Typically, the cross-talk increases the error per 
observation of the third-star from the expected level (for a sharp-line star) of 1.2 — 1.3 km s~^, to 
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the level of 1.5 — 2.5 km s . Except for noting the presence of the cross-talk, we are not in position 

to study it more extensively given the different type of the binary-phase dependence in each case. 
Our hope is that the cross-talk will average out in the velocities of the third component, although 
the final proof will come only through external comparisons. The stress has been always on the 
quality of the binary solutions, perhaps at the expense of the quality of the radial-velocity data for 
the third components. 

We measure the radial-velocities for the binary components, and - if necessary, of a spectro- 
scopic companion - but do not estimate the accuracy of the radial-velocity measurements at this 
stage. In principle, it is possible to establish a relation between the signal-to-noise (S/N) in the 
spectrum and in the broadening function (Rucinski et al. 1993), but further propagation of the 
errors into the velocity errors is more complex and depends on many factors. While such a rela- 
tion would be definitely needed for full modeling of the BF's, we feel that complexity of the error 
analysis is not warranted in our case. Thus, we do not determine the radial-velocity errors from 
the individual BF's, but evaluate them externally later from the orbital velocity solutions. Such 
estimates may be perhaps overly pessimistic, as they incorporate systematic deviations from the 
orbital motion models. Most importantly, however, the random errors are not the limiting factor 
in our results; the real difficulty is in evaluation of systematic uncertainties. Wc address this issue 
in Section 9, after describing the orbital solutions (Section 7) and the externally-evaluated random 
errors (Section 8). 



7. ORBITAL SOLUTIONS 

Step #4 of the reductions is the determination of the radial-velocity orbit using individual 
velocities of both components at all observed orbital phases. Currently, we do so by measuring 
individual velocities of components, although a more global approach involving modeling of broad- 
ening functions would be definitely much preferable. The broadening functions have the potential 
of providing much more information than just simple velocity centroids so that orbital solutions 
could be carried to a much higher level of sophistication than in this series of papers. Such use 
of the BF's was described in Rucinski (1992); Lu & Rucinski (1993); Rucinski et al. (1993), where 
modeling of the BF shape was advocated. Full modeling of this type requires knowledge of the 
orbital inclination, which is usually not available, and involves a simultancoTis determination of 
the radial- velocity span, Ki + K2, the mass ratio, q and the degree-of-contact parameter, /. The 
complexity of such a global approach is the main reason why we continue to use single velocities to 
characterize motions of stellar components, but we do recognize limitations of this approach which 
may generate systematic uncertainties in the final results; this is discussed in Section 9. We should 
add, that originally, this program was intended to provide the Vq values from a small number of 
radial-velocity measurements, to relate to the then newly available Hipparcos tangential velocities. 
However with time, our program acquired its current significance as the main contributor of ra- 
dial velocity orbits for short-period binaries (this circumstance taking place partly "by default". 
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through a surprising lack of similar programs at other observatories). Thus, we continue to use the 
Gaussian fits, but recognize that all our spectra and the broadening functions may be used for a 
much more extensive modeling. 

All short-period binary systems observed by us so far have circular orbits resulting in sine- 
curve variations of orbital velocities. The only exception that we had to consider is the third 
star in the system of HT Vir (Paper IV) on an eccentric orbit; for this case we used the model 
of Morbey (1975). Because eclipse effects of rotationally broadened lines change line shapes and 
produce undesirable radial-velocity shifts, we eliminate observations close to orbital conjunctions, 
usually within the phase ranges 0.85 - 0.15 and 0.35 - 0.65. 

The orbital solutions are obtained iteratively. First, we use the linear model of two sine 
curves and one constant value, with an assumed moment of the primary eclipse Tq. Thus, for k 
observations, we simultaneously fit by least squares 2 x A; equations of the type: 

Vi{(j)i) = Vo- Ki sin cl>i + (5) 
V2{(t>i) = Vo- +K2 sin (1)1 (6) 
I = 0,...,k-l 

is the orbital phase, (fn = [ti — Toj/P. Similarly to Tq, the period P is usually taken from literature 
sources and is fixed; only in a few cases we attempted to improve its value. The equations can 
be weighted at this point when observations are of different quality. The weighting schemes are 
discussed in descriptions of stellar systems in individual papers and are given in the tables with 
radial velocities. Note the sign convention in the equations, which implies that we usually start with 
an assumption that the primary, more massive component (star 1) is eclipsed at the photometric 
primary minimum. In other words, we assume that for a contact system the configuration is of an 
A-type contact binary. We identify the W-typc systems when this assumption is not valid. 

The resulting Vq, Ki, K2 are the first approximations of the orbital parameters. The next 
step in the iterative solution consists of the application of the linearized versions of Equations (5) 
- (6) for AVb, AKi, AK2 and ATq. We always first use any available literature value for Tq and 
then improve it by solving the linearized equations until all corrections A no longer change. It 
is at this stage that we determine random-error uncertainties of the orbital parameters and the 
radial-velocity errors per observation. 



8. MEAN STANDARD ERRORS 

Least-squares solutions of the linearized Equations (5) - (6) can provide mean standard errors 
of the orbital parameters Vq, Ki, K2 and Tq. We do not use such errors because they usually 
underestimate the random error uncertainties. Instead, we use the "bootstrap sampling" technique 
which involves multiple (thousands of times) re-sampling of the data with possible repetitions, 
with subsequent solutions of all such datasets. By forming statistics of the spread in the resulting 
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parameters and by determining the inner 67 percent distribution ranges, we estimate equivalents of 
the mean standard errors. They are sometimes close to the linear estimates, but usually are larger. 
In any case, we consider them more realistic as they include inter-parameter correlations. 

We have sufficient amount of information from all our orbital solutions to analyze sizes and 
distributions of our random errors. For that purpose, we used all the available solutions, eliminating 
three systems observed at the Dominion Astrophysical Observatory, as reported in Paper I, and 
adding W Crv described separately (Rucinski & Lu 2000), altogether 58 orbital solutions. 

The statistics of mean standard errors per single observation, Cj, for the primary (i = 1) and 
the secondary (i = 2) components, is shown in the first panel of Figure 3. The median values of 
the errors are (ei) = 5.48 km s"-*^ and (€2) = 11.50 km s"'^. The corresponding distributions for the 
errors of the radial-velocity amplitudes, a{Ki), are shown in the upper right panel. The median 
values are {a{Ki)) = 1.11 km s~^ and {a{K2)) = 1.96 km s~^. The center-of-mass velocities 
Vo arc better established than Ki because two stars contribute in each solution to one number. 
The median value of these errors is ((t(Vo)) = 1.07 km s~^. Finally, the distribution of the mean 
standard errors of the initial epoch, a{To), is shown in the last panel of Figure 3. The median value 
for this error is, {a{To)) = 0.0011 days (about 1.5 minute). 

The mean standard errors of the orbital parameters are correlated. The most interesting 
correlations are shown in Figure 4. The two upper panels show that the mean error of the center- 
of-mass velocity, o"(Vb), which appears to be a convenient measure of the quality of the orbital data. 
It depends on the brightness of the system and on the orbital period. It is confined within < 1.5 
km s^^ for Vmax < 8.5, but increases to slightly over 2 km s^-*^ for Vmax > 10 (the left upper panel). 
The scatter in a(Vb) increases for short-period systems (the right upper panel of Figure 4), but this 
may be due to the fact that most of our targets had periods within 0.3 - 0.6 days, in a range where 
a genuine frequency maximum exists in the volume-limited samples of contact binaries (Rucinski 
1998). The systems with longer periods, P > 0.8 days, tend to show small errors, but these are 
exactly those binaries which had been overlooked before among bright stars and have been easy 
targets for our program. The error o"(Vo) correlates tightly with a{Ki) and with e.j, as shown in the 
four lower panels of Figure 4. A particularly close correlation with the slope close to unity exists 
between a{Ki) and cr{Vo) (middle left panel). 

Binaries observed in spectroscopic triple systems show slightly larger random errors than when 
isolated, typically by less than 1 km s~^ in the errors of Vq, Ki and K2. Such binaries are shown 
by open symbols in all plots in Figure 4. 

9. SYSTEMATIC UNCERTAINTIES 

It is difficult to evaluate systematic uncertainties of our results. The systematic errors depend 
in a complex way on the orbital parameters and couple with the random errors. The main source 
of the systematic errors is the measurement of radial velocities from the broadening functions. We 
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approximate the center-of-mass positions with the hght centroids and measure the centroids by 
fitting Gaussians. The latter assumption, that the radial velocities of the light centroids coincide 
with radial velocities of the mass centers, is - in general - not fulfilled by distorted components in 
close binary systems and is particularly dangerous for contact binaries where the peaks in the BF's 
are not symmetric, with steeper outer parts and more gently sloping inner parts. Direct modeling 
of the BF's would avoid this systematic error (see the end of this section). 

Some insight into systematic uncertainties involving the Gaussian approximation of the peaks 
in the broadening functions can be obtained by applying Gaussians of various widths and evaluating 
systematic shifts in the results. We will consider here, as a case study, of a typical, 9 magnitude, 
A-type contact system GM Dra from the most recent Paper VI. 

Let us first consider one broadening function for the orbital phase 0.283 of GM Dra (the lowest 
panel of Figure 5). For this particular broadening function, we would normally select the best-fitting 
Gaussians to have the width parameters 031 = 120 km s^-*^ and 032 = 80 km s~-^ (see Eq. (4)). 
However, as an experiment, we considered widths between the estimated narrowest and widest 
acceptable values of a^f 031 = 100 — 140 km and 033 = 60 — 100 km s~^. The extreme cases are 
shown by dotted and broken lines in Figure 5. For the full ranges of the widths, the change in the 
measured velocity of the primary component is from —28.86 km to —31.27 km s~^, while the 
change for the secondary component is from +261.44 km s~^ to +262.49 km s~^. Thus, systematic 
errors in radial velocities appear to be at a level of 1.5 to 2.5 km s^^, with larger velocities (in the 
absolute sense) associated with larger assumed widths of the fitting Gaussians. 

Analysis of the type presented above can be done for all available broadening functions, at all 
orbital phases. The four upper panels of Figure 5 show the shifts in the measured centroids for 
all available observations of GM Dra obtained around the orbital quadratures, within the orbital 
phase ranges 0.15 - 0.35 and 0.65 - 0.85, as marked in the figure. The Gaussian widths 031 were 
incremented in equal steps and for each assumed width a full radial velocity determination was 
performed. As we can see in the figure, the systematic effects are clearly present, especially for 
the secondary (less-massive) component. The shifts are typically at the level below 2 km s^^ for 
the primary component, but the shifts of the order of 5 - 7 km s~^ are not uncommon for the 
secondary component. The shifts depend on side of the binary system (or the sign of the radial 
velocity) observed at a given orbital quadrature. The overall tendency appears to be that the wider 
Gaussian width 03^ result in velocities further away from the center-of-mass velocity, i.e. such ones 
which should lead to systematically larger values of the orbital amplitudes, Ki. This is confirmed 
by the actual determinations of the radial velocity orbits for the extreme values of [031, 032] pairs, 
selected to deviate from the optimal values of 120 and 80 km s~^ by ±20 km s~^. The systematic 
changes for the particular case of GM Dra strongly depend on the parameter considered. While 
the changes in Vq are within +0.09 and —0.16 km s~^, those in the amplitudes are larger: —0.34 
and +0.15 km for Ki and as much as —4.30 and +5.97 km for K2. While the ranges of 
the Gaussian widths 0,3,4 were intentionally exaggerated in the experiment, to estimate the largest 
systematic deviations, we clearly see that systematic effects may set an important limitation on 
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our results. For comparison, we note that the random errors of the orbital parameters of GM Dra 
are a{Vo) = 1.52 km s'^, a{Ki) = 1.75 km s'^ and a{K2) = 2.50 km s'^ (Paper VI). Thus, for 
this particular binary, the systematic uncertainty appears to be larger than the random error only 
for K2, but then it is even two times larger. 

Optimally, the systematic effects resulting from the use of different widths in the Gaussian fits 
should be evaluated for each binary through a process similar to that applied to GM Dra. However, 
we feel that it is impractical to perform similar analyses for all systems in this program. Besides, 
we know that application of the Gaussian fits is - in any case - a crude approximation and that 
the best approach would be to model the broadening functions as it was done in Rucinski (1992); 
Lu & Rucinski (1993); Rucinski et al. (1993). Full BF modeling would permit inclusion of more 
spectra than we utilize now because, currently, we measure for radial velocities only those BF's 
which show a clear splitting of the spectral signatures. By addition of these spectra we would 
increase the available material by about 20 - 30 percent, which would only slightly reduce random 
errors and thus produce a very modest improvement in accuracy. Much more important would be 
a reduction or entire elimination of the systematic errors, which may reach levels of 5 - 7 km s~^. 
For most binaries of this program, this would typically correspond to about 2-3 percent error 
in Ki, but in some extreme cases of small semi-amplitudes, the errors may reach 10 - 15 percent. 
While the approach involving combined radial velocity and light-curve modeling would avoid the 
main systematic effects, it would require a considerable organizational and computational effort, 
introducing large delays in our mostly observational program. Since our radial velocity observations 
are — for most systems - the first and the only ones, we decided to accept the level of systematic 
errors generated by the use of the measuring Gaussians and make our solutions generally available, 
keeping in mind their systematic uncertainties which must be taken into account when considering 
the overall accuracy of our program. 

10. CONCLUSIONS AND PLANS 

The ongoing survey of close binary systems with periods shorter than one day, currently 
conducted at the David Dunlap Observatory, has resulted in a consistent set of radial velocity 
orbits for sixty previously unobserved binaries to approximately 11th magnitude. While, at the 
start, the survey concentrated on systems which simply had not been studied before (for various 
reasons, but mostly because of inadequate instrumentation and data-analysis tools some half a 
century ago, when this field was very active), the photometric discoveries of the Hipparcos satellite 
are now dominating in numbers. There was only one Hipparcos system among the first twenty orbits 
(Papers I and H), 9 such systems among the next twenty orbits (Papers HI and IV) and 15 such 
systems among the most recent twenty orbits (Papers V and VI) . About 50 known, photometrically- 
discovered binaries still remain to be observed and analyzed and new ones are constantly added 
to catalogs, some of them quite bright. Regrettably, apparently there is no similar survey for the 
southern hemisphere. 
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Our survey is quasi-random in the sense that we observe all short-period (P < 1 day), bright, 
previously unobserved binaries. With such criteria, the contact binaries absolutely dominate in 
numbers. Among the 60 systems described in the previous six papers, only 8 were not contact 
systems. This is partially due to strong selection effects against detection of detached binaries, 
but mostly due to the very high frequency of contact binary systems in the old-disk population, 
particularly in the period range 0.3 to 0.5 days, but with a tail extending beyond one day, to about 
1.3 - 1.5 days. The high frequency of incidence is strongly manifested in the volume- limited OGLE 
sample and in open clusters (Rucinski 1998). Because our survey is magnitude limited, we tend to 
include many brighter systems from the tail of the distribution between 0.5 day and our current 
upper limit at one day. Otherwise, we do not discriminate among binary systems in any other 
way. In particular, the random character of the survey has resulted in discoveries of the largest 
{q = 0.97, V753 Mon; Paper III) and the smallest (g = 0.066, SX Crv; Paper V) known mass ratios 
among contact binaries. 

The DDO survey is characterized by moderate random errors of about 1-2 km s"^ for 
the orbital parameters, Vq, Ki and K2, and - upon completion - can serve as a useful database 
of parameters of very close binary systems. We are aware, however, that our final parameters 
contain systematic uncertainties resulting from our radial- velocity measurement techniques. While 
the use of the broadening functions permitted us to analyze close binaries in several multiple, 
visual/spectroscopic systems providing data which were too "difficult" before, our extraction of 
individual radial velocities from the broadening functions, through Gaussian fitting, is a disputable 
approach for contact binary systems. Because the line-broadening for such systems is very strong, 
comparable with orbital velocities of hundreds of km s^^, and - in fact - somewhat asymmetric, 
our measuring technique may lead to systematic errors reaching levels of 5 - 7 km s^^ or even 
more. Paradoxically, through the use of the broadening functions in place of the cross-correlation 
functions, we have uncovered real physical reasons why the Gaussian approximation is only barely 
appropriate. The correct approach avoiding the systematic errors would be to model the broadening 
functions and determine the radial velocities in terms of the mass ratio, q, and the scaling factor, 
(K1+K2), with the shift, Vq. The models would require independent input from parallel solution of 
light curves, providing the orbital inclination angle, i, as well as the degree-of-contact, /. Currently, 
most of the program targets have not had their light-curves solved, and even if some attempts have 
been made, we would not trust them for the following simple reason: We have seen so many cases 
of the spectroscopic mass ratio different from the previous photometric mass-ratio determinations, 
Qsp Qphot, that we feel very strongly that the values of qphot are usually not properly constrained 
and may be plainly wrong^, leading to entirely incorrect combinations of orbital parameters. 

We envisage that the results of this survey will provide just a first stage of an iterative process. 



^Totally eclipsing systems are an exception, as pointed by Mochnacki & Doughty (1972a, b), but then chances of 
total eclipses depend on the mass-ratio itself (a wider range of inclinations for small values of q), producing a very 
complex bias in the uncertainties of qphot- 
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In future, our spectroscopic values of mass ratio, g^p, should permit solution of light curves which 

were previously unsolvable because to the poorly constrained mass ratios. The derived information 
on (i, f) pairs would permit, in turn, a re-discussion of the broadening functions and determination 
of the final orbital parameters, free of systematic uncertainties. 

Concerning the instrumental developments at the DDO: Soon, we plan to start using a new 
CCD system based on a much more sensitive detector. While the analysis of the data should remain 
the same as described above, we may have to select the targets more discriminately. In particular, 
it may turn out impractical to observe all binaries with periods shorter than one day down to to 
the expected limiting magnitude of about 12.5 magnitude. Indeed, from the point of astrophysical 
usefulness, it would be advantageous to reduce the deficit of the intrinsically faint contact systems 
among spectroscopically studied binaries of the magnitude-limited sample, by attempting to form 
a volume-limited sample through giving preference to very short-period systems. 

While many persons have participated in this program and have either co-authored the previous 
papers or their contributions have been acknowledged there, special thanks are due to Dr. Hilmar 
Duerbeck who contributed to setting the goals of the program in its early stages when it was 
concerned mostly with the center-of-mass velocities for contact binaries for a planned spatial- 
velocity investigation. 

The author would like to thank Stefan Mochnacki and Mel Blake for reading and commenting 
on an early version of the paper. Thanks are also due to the anonymous referee for two very careful 
and constructive reviews. 
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Captions to figures: 

Fig. 1. — The left panel shows spectra of the sharp-lined template star, HD 89449 (upper spectrum), 
and of the close binary star KR Com (lower spectrum, shifted down by 0.5 in the observed flux). 
The spectral types of the stars are F6IV and GOIV, respectively. The right upper panel shows 
the broadening function (BF) obtained by our linear de-convolution using the two spectra in the 
left panel. When measuring radial-velocities of a binary, we initially fit the whole triple feature 
by Gaussians, then subtract the sharp-line component and repeat the determination for the close 
binary. The three components of the BF are shown by thin lines. The lower part of the right panel 
shows the cross-correlation function (CCF, thick line), in comparison with the BF (thin line), both 
obtained from the same spectra at left. The CCF has much lower resolution than the BF, but also 
shows negative excursions in the zero (baseline) level. While Gaussians may be a reasonable tool 
for measurement of radial velocities from the CCF's, the BF's are much better defined; note the 
much steeper outer ends of the BF relative to the Gaussians. We discuss extensively the systematic 
uncertainties of this type in Section 9. 

Fig. 2. — This figure shows an experiment in data processing. In the left panel, a high-resolution 
spectrum, rebinned to equal velocity steps of 0.88 km s~^, is shown without any additional broad- 
ening (dotted line) and with rotational broadening of V sin i = 88 km s"^. The CCF for the two 
spectra is shown in the right panel. Notice the strong positive fringes outside the main correlation 
peak as well as the shift of the quasi-baseline well below the expected zero level; the actual broad- 
ening function has been shifted down by —0.1 units to visualize the most likely placement of the 
local baseline in the vicinity of the correlation peak. 

Fig. 3. — Distributions of mean standard errors for program binaries. The histograms give the 
distributions of the error per observation (for each component), e^, and of the errors of orbital 
parameters a{Vo), (t{Ki) and a{K2), all expressed in km s~^. The bin sizes are given in the y-axis 
labels. In the two upper panels, full-line histograms are for the primary components (ei and a{Ki)) 
while the dotted histograms are for the secondary components (e2 and a{K2)). The last panel gives 
the distributions of mean standard errors for the initial epoch Tq in units of 0.0001 days. 

Fig. 4. — Correlations between various mean standard errors, as given in the axis labels. The two 
upper panels show cr(Vb) as a function of Vmax and the orbital period, P. cr(Vb) is a convenient 
measure of the solution quality and correlates tightly with a{Ki), as shown in the middle left panel. 
The middle right panel shows the correlation between (7(Ki) and a{K2). This correlation is not 
perfect because of the very large range of mass ratios observed among binaries of this program. The 
two lowest panels show the errors per observation, e; the left panel shows the correlation between 
ei and <y{Ki) while the right panel shows the correlation between e\ and €2- In all panels binaries 
analyzed through subtraction of the third-component signatures from the broadening functions are 
marked by open circles. All quantities are expressed in km s~^. 

Fig. 5. — The four upper panels show shifts in the measured velocities for the primary (left panels) 
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and secondary (right panels) components of GM Dra (Paper VI) versus the Gaussian-width pa- 
rameter a^i (see text). Each hne is for one broadening function at one orbital phase within ranges 
around the two orbital quadratures, as marked in the panels. The lowest single panel shows one 
broadening function of GM Dra at the orbital phase 0.283, approximated by the Gaussians with 
the width parameters [a^i, 032] considered most extreme for this case: [100, 60] km s~^ (dotted 
line) and [140, 100] km s~^ (broken line). 
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